#!/usr/bin/env python

" modele d'ising bidimentionnel "

from pygame import *
from math import *
from random import *
import sys

def notice() :
	" notice du programme "

	print ""
	print " ##########################################################"
	print " ##            Modele d'Ising bidimensionnel             ##"
	print " ##                 Transition de phase :                ##"
	print " ##           ferromagnetisme - paramagnetisme           ##"
	print " ##########################################################"
	print ""
	print " Cliquez sur les grilles pour initialiser le materiau dans l'une"
	print " des configurations suivantes :"
	print "  -en damier : energie maximale, tous les spins voisins sont retournes"
	print "  -aleatoire : energie nulle, magnetisation nulle"
	print "  -spins paralleles : energie minimale, magnetisation maximale"
	print ""
	print " Cliquez sur les boutons bleux pour changer le nombre d'etats accessibles pour chaque composante",
	print "ainsi que le nombre de degres de liberte pour chaque spin"
	print ""
	print " Cliquez sur le thermometre pour faire varier la temperature."
	print ""
	print " Cliquez sur le bouton Stop, pour arreter ou redemarrer l'animation"
	print ""
	print " Utilisez le bouton 'H' pour faire varier le champ magnetique exterieur"


# constantes du programme

	# nombre de spins modifies par refresh			 			
STEP  = 500

	# temperature maximum		
TEMP_MAX = 10


# parametres d'affichage

	# rayon des spins, bordure
R, DECAL, GRID = 5, 20,  [70, 70]

	# grilles
iX, iY, iW, iH = GRID[0] * R + 2 * DECAL, DECAL, 7 * R, 7 * R

	# spins
sX, sY  = GRID[0] * R + 2 * DECAL, 2 * DECAL + 7 * R, 
sY2 = sY + 2 * DECAL

	# thermometre
tX, tY, tW, tH = GRID[0] * R + 2 * DECAL + 5, sY + 4 * DECAL, 10, 140

	# effets de bords
eX, eY, eW, eH = GRID[0] * R + 3 * DECAL + 5 + tW, sY + 6 * DECAL, 15, 15

	# stop and play
bX, bY, bW, bH = GRID[0] * R + 2 * DECAL, tY + tH + 2 * DECAL , 20, 20

	# champ H
cX, cY = bX + bW + 2 * DECAL + 60, bY - 5	

	# graphes
g1X, g1Y, g1W, g1H = DECAL, max(bY + 2 * DECAL, R * GRID[1] + 2 * DECAL), 200, 100
g2X, g2Y, g2W, g2H = g1X + g1W + DECAL, g1Y, g1W, g1H

	# taille de la fenetre
WIN_W = max(GRID[0] * R + 3 * DECAL + 6 * R + 6 * iW, g1W + g2W + 5 * DECAL)
WIN_H = g1Y + g1H + DECAL


class Graph :
	" genere un graphique anime "

	def __init__ (self, posx, posy, Width, Height, legende):
		
		# position du point sur le graph
		self.cursor = 0

		# largeur du graph
		self.posx = posx
		self.posy = posy
		self.Width = Width
		self.Height = Height
		self.curs_height = Height
		self.legende = legende

	def trace(self, screen, y, value_max):
		" met a jour le graphique "
		
		# calcule les coordonnees
		self.coord = [self.posx + self.cursor, 
									self.posy + self.curs_height / 2 - y * self.curs_height / (2 * value_max)]
	
		if self.cursor == 0:
			# affiche le fond du graphique
			draw.rect(screen, [100, 100, 100], 
										(self.posx, self.posy , self.Width, self.Height + 1), 0)
										
			# trace l'abscisse y = 0					
			draw.line(screen, [0, 0, 0], 
										(self.posx, self.posy + self.Height/2), 
										(self.posx + self.Width, self.posy + self.Height/2))
			self.last = self.coord
			
			# affiche la legende
			f = font.SysFont("Times New Roman", 13)
			s = f.render(self.legende, False, (0, 0, 0))
			screen.blit(s, (self.posx+10, self.posy+10))
			
	 
	 	# trace le segment qui joint l'ancien point au nouveau
		draw.line(screen, [255, 255, 255], self.last, self.coord)
		
		# memorise le point actuel
		self.last = self.coord
		
		# incremente le curseur
		self.cursor = (self.cursor + 1) % self.Width
		
												
class IsingAnim:
	" animation "
	
	def __init__ (self, spin):
		" initialisation de l'animation "
		
		# creation de la fenetre
		init()		
		self.screen = display.set_mode([WIN_W, WIN_H])
		
		# initialisation de la grille, des boutons et du thermometre
		self.init_grid(spin)
		
		# creation des boutons d'initialisation
		font.init()
		
		# grille en damier
		bit, bit2 = 1, 1
		for i in xrange(7):
			bit2 *= -1
			bit = bit2	
			for j in xrange(7):
				bit *= -1
				if bit > 0: col = [150, 200, 0]
				else: 		col = [0, 0, 200]
				draw.rect(self.screen, col, (5*i+iX, 5*j+iY, 4, 4), 0)
		
		# grille aleatoire
		for i in xrange(7):
			for j in xrange(7):
				col = choice([[150, 200, 0], [0, 0, 200]])
				draw.rect(self.screen, col, (5*i+iX + 10 * 5, 5*j+iY, 4, 4), 0)
		
		# grille verte		
		col = [150, 200, 0]
		for i in xrange(7):
			for j in xrange(7):
				draw.rect(self.screen, col, (5*i+iX + 20 * 5, 5*j+iY, 4, 4), 0)
		
		# boutons et thermometre
		self.buttons_up ('play', 1)
		self.buttons_up ('field', None)
		self.thermo_up(tY+tH-tH/2-2, TEMP_MAX/2)

		
		# creation des graphiques
		self.graph_E = Graph(g1X, g1Y, g1W, g1H, "Energie moyenne")
		self.graph_M = Graph(g2X, g2Y, g2W, g2H, "Magnetisation moyenne")		
	
	
	def init_grid (self, spin):
		" affiche la grille de spins "
		
		for i in xrange(GRID[0]):
			for j in xrange(GRID[1]):
				self.change_spin(i, j, spin[i][j])
		self.refresh()
	
	
	def buttons_up (self, but, state, state2 = 0):
		" mise a jour des boutons	"
		
		if but == 'play':
			col1, col2 = [0, 0, 100], [100, 0, 0]
			
			draw.rect(self.screen, col2, (bX, bY, bW, bH), 0)
			if state == 0: draw.lines(self.screen, col1, 1,
				((bX+5, bY+5), (bX+bW-5, bY+5), (bX+bW-5, bY+bH-5), (bX+5, bY+bH-5)))
			else: draw.lines(self.screen, col1, 1,
				((bX+5, bY+5), (bX+15, bY+10), (bX+5, bY+15)))
				
		elif but == 'spin':
			if state == 0:	sp = 'X'
			else:						sp = str(state)
			draw.rect(self.screen, [0, 0, 0], (sX, sY,190, bH) ,0)
			draw.rect(self.screen, [0, 0, 255], (sX, sY,bW, bH), 0)
			
			f = font.SysFont("Times New Roman", 15)
			s = f.render(sp, True, (0, 0, 0))
			self.screen.blit(s, (sX+5, sY-1))
			
			s = f.render('etats par composantes', True, (0, 0, 255))
			self.screen.blit(s, (sX+bW+10, sY-1))

		elif but == 'tmc':
			draw.rect(self.screen, [0, 0, 0], (tX + 30, tY +80, 300, 20) ,0)
			
			f = font.SysFont("Times New Roman", 15)
			s = f.render("temps de monte carlo : %d" % (state), True, (0, 0, 255))
			self.screen.blit(s, (tX + 30, tY+80))

		elif but == 'bord':
			draw.rect(self.screen, (0, 0, 0), (eX+1, eY+1, 200, 30), 0)
			f = font.SysFont("Times New Roman", 15)
			s = f.render('Effets de bords', True, (0, 0, 255))
			self.screen.blit(s, (eX+eW+DECAL, eY))
			draw.rect(self.screen, (0, 0, 255), (eX, eY, eW, eH), 0)
			draw.rect(self.screen, (0, 0, 0), (eX+1, eY+1, eW-2, eH-2), 0)

			if state:
				draw.lines(self.screen, (0, 0, 255), 1, 
										((eX+1, eY+1), (eX+eW-1, eY+eH-1)))
				draw.lines(self.screen, (0, 0, 255), 1, 
										((eX+1, eY+eH-1), (eX+eW-1, eY+1)))
		
		elif but == 'liberty':
			sp = str(state)
			
			draw.rect(self.screen, [0, 0, 0], (sX, sY2,140, bH) ,0)
			draw.rect(self.screen, [0, 0, 255], (sX, sY2,bW, bH), 0)
			
			f = font.SysFont("Times New Roman", 15)
			s = f.render(sp, True, (0, 0, 0))
			self.screen.blit(s, (sX+5, sY2-1))
			
			s = f.render('degres de liberte', True, (0, 0, 255))
			self.screen.blit(s, (sX+bW+10, sY2-1))
			
		elif but == 'field':
			draw.circle(self.screen, [0, 0, 0], (cX, cY), 40)
			draw.rect(self.screen, (0, 0, 0), (bX+bW+DECAL, bY-15, 40, 30), 0)
			f = font.SysFont("Times New Roman", 15)
			s = f.render('H', True, (0, 0, 255))
			self.screen.blit(s, (bX+bW+DECAL, bY))
			draw.rect(self.screen, (0, 0, 255), (bX+bW+DECAL+16, bY+5, 15, 15), 0)
			draw.rect(self.screen, (0, 0, 0), (bX+bW+DECAL+17, bY+6, 13, 13), 0)
			
			if state:
				draw.lines(self.screen, (0, 0, 255), 1, 
										((bX+bW+DECAL+16, bY+5), (bX+bW+DECAL+29, bY+18)))
				draw.lines(self.screen, (0, 0, 255), 1, 
										((bX+bW+DECAL+17, bY+19), (bX+bW+DECAL+30, bY+5)))
				draw.circle(self.screen, [0, 0, 255], (cX, cY), 40)
				draw.circle(self.screen, [0, 0, 0], (cX, cY), 39)
				draw.circle(self.screen, [0, 0, 255], (cX, cY), 0)		
				
				if state2 :
					 R, V, B = 0, 0, 0
			 		 R = (state2[0] + 4) * 255 / 8
			 		 V = (state2[1] + 4) * 255 / 8
			 		 B = 127
			 		 draw.lines(self.screen,[R, V, B], 1, ((cX, cY), state))
			
		self.refresh()

		
	def thermo_up (self, y, temp):
	
		" change l'etat du thermometre "
		draw.circle(self.screen, [150, 0, 0], (tX+tW/2, tY+tH+12), 14)
		draw.rect(self.screen, [0, 0, 100], (tX-1, tY, tW+1, tH), 0)
		draw.rect(self.screen, [150, 0, 0],	(tX, tH + tY, tW, y - tY - tH), 0)
									
		draw.rect(self.screen, [0, 0, 0], (tX + 30, tY, 100, 20), 0)
		f = font.SysFont("Times New Roman", 15)
		s = f.render("Temp: %.2f" % (temp), True, (0, 0, 255))
		self.screen.blit(s, (tX + 30, tY))
		
		self.refresh ()

		
	def change_spin (self, i, j, s, grille = 1):
		" change le spin (i, j) de la grille "
		
		# calcule la couleur en fonction du spin
		col = [(s[0]+1)*255/2,(s[1]+1)*255/2,(s[2]+1)*255/2]
		
		# dessine un carre de rayon r avec la couleur voulue
		draw.rect(self.screen, col, (R * i + DECAL, R * j + DECAL, R-grille, R-grille), 0)

		
	def refresh (self):
		" refresh de l'animation "
		
		display.update()

		
		

class IsingModel:
	" modele d'Ising "
	

	def __init__ (self, style, display = 1, nb_etat = 2, liberty = 1, H = [0,0,0], J = -1, grille = [70,70]):
		" initialisation "
		
		global GRID, NB_SPIN
		
		# taille de la grille	
		GRID = grille

		# nombre de spins
		NB_SPIN = GRID[0]*GRID[1]

		# initialise la temperature
		self.temp = TEMP_MAX / 2

		# champ magnetique exterieur
		self.H = H
		
		# constante de couplage						 	
		self.J = J
				
		# 2 ou X etats par spins?
		self.nb_etat = nb_etat
		
		# nombre de degres de liberte (1, 2 ou 3)
		self.liberty = liberty
	
		# l'animation est inactive		
		self.stop = 1
		
		# genere la grille
		self.init_grid(style)
		
		# redirection des resultats sur l'ecran?
		self.display = display

		# effets de bord
		self.bord = 0

		# Si oui
		if self.display :
			
			# on cree l'animation
			self.anim = IsingAnim(self.spin)

			# on actualise les boutons
			self.anim.buttons_up ('spin', self.nb_etat)
			self.anim.buttons_up ('liberty', self.liberty)
			self.anim.buttons_up ('bord', self.bord)

			# on affiche le mode d'emploi
			notice()

			 
	def get_new_spin(self, s = [0,0,0] ):
		" renvoit un nouvel etat "

		# on tire au hasard un etat different parmi les etats disponibles parmi 2 etats
		if self.nb_etat == 2:
			
			# 1 seul degre de liberte
			if self.liberty == 1: 
				while 1 :
					s2 = choice([[1,0,0],[-1,0,0]])
					if s2 != s : break

			# 2 degres de liberte
			elif self.liberty == 2:
				 while 1 :
					 s2 = choice([[1,0,0],[-1,0,0],[0,1,0],[0,-1,0]])
					 if s2 != s : break
		
			# 3 degres de liberte
			elif self.liberty == 3:
				 while 1 :
				 		s2 = choice([[1,0,0],[-1,0,0],[0,1,0],[0,-1,0],[0,0,1],[0,0,-1]])
				 		if s2 != s : break
		
		# parmi un etat d'orientation quelconque de norme <= 1
		else :
			while 1:
				norme = 0
				var = [0 for i in xrange(3)]
				for i in xrange(self.liberty):
					var[i] = uniform(-1, 1)
					norme += var[i]**2
				if norme <= 1 : break
			s2 = [var[0],var[1],var[2]]
		
		return s2

		
		
	def init_E_M(self) :
		" on initialise l'energie et la magnetisation totale du systeme	"
		
		# magnetisation totale
		self.M_total = .0
		
		# energie totale
		self.E_total = .0


		# parcours de la grille
		for i in xrange(GRID[0]) :
			for j in xrange(GRID[1]) :
			
				# calcule la magnetisation totale
				self.M_total += self.spin[i][j][0] + self.spin[i][j][1] + self.spin[i][j][2]

				# calcule l'energie totale
				for n in xrange(3) :
					self.E_total += self.spin[i][j][n] * (self.J * (self.spin[(i-1)%GRID[0]][j][n] + self.spin[i][(j-1)%GRID[1]][n]) - self.H[n])

		
	def init_grid (self, style='hot'):
		" initialisation de la grille "
		
		# initialise la grille de spins
		self.spin = [ [[0,0,0] for j in xrange(GRID[1])] 
			 						 for i in xrange(GRID[0]) ]
		
		# Distribution aleatoire
		if style == 'alea' :	
			for i in xrange(GRID[0]):
				for j in xrange(GRID[1]):
						self.spin[i][j] = self.get_new_spin()

			 
		# Distribution en damier, energie maximisee 
		elif style == 'hot' :
			s, s2 = 1, 1
			for i in xrange(GRID[0]):
				s2 *= -1
				s = s2
				for j in xrange(GRID[1]):
					s *= -1
					self.spin[i][j] = [s,0,0]

	 
		# Spins tous paralleles
		elif style == 'cold' :
			for i in xrange(GRID[0]):
				for j in xrange(GRID[1]):
						self.spin[i][j] = [-1,0,0]

		# initialise l'energie et la magnetisation totale							 
		self.init_E_M()
									 
 

	
	def event_handler (self, events):
		" gestion des evenements "
		
		for ev in events:
			if ev.type == QUIT:
				sys.exit()
			if ev.type == MOUSEBUTTONDOWN:
				
				# thermometre
				if tX <= ev.pos[0] <= tX + tW and tY + 2 <= ev.pos[1] <= tY + tH + 24:
						
					# changement de la temperature (echelle: 0.01 - 10)
					if ev.pos[1] > tY + tH:
						self.temp = .01
					else:
						self.temp = float(tY + tH / 2 - ev.pos[1] + 71) * TEMP_MAX / tH

					self.anim.thermo_up(ev.pos[1], self.temp)
				
				# effets de bord
				elif eX <= ev.pos[0] <=  eX + eW and eY <= ev.pos[1] <= eY + eH:
					self.bord = (self.bord + 1) % 2
					self.anim.buttons_up('bord', self.bord)

				# selectionne le champ magnetique exterieur
				elif cY-40 <= ev.pos[1] <= cY+40 and cX-40 <= ev.pos[0] <= cX+40:
				
					if sqrt((cY-ev.pos[1])**2 + (cX-ev.pos[0])**2) > 39: break
			
					if not (self.H[2]): break
			
					self.H[0] = float(ev.pos[0] - cX) / 10
					self.H[1] = float(cY - ev.pos[1]) / 10
					
					# recalcule l'energie totale
					self.init_E_M()

					self.anim.buttons_up('field', (ev.pos), self.H)

				# selectionne le nombre d'etats
				elif sY <= ev.pos[1] <= sY+bH and sX <= ev.pos[0] <= sX+bW:
					self.nb_etat = (self.nb_etat + 2) % 4
					self.anim.buttons_up('spin', self.nb_etat)
				
				# selectionne le nombre de degres de liberte
				elif sY2 <= ev.pos[1] <= sY2+bH and sX <= ev.pos[0] <= sX+bW:
					self.liberty = (self.liberty + 1) % 3
					if self.liberty == 0 : self.liberty = 3
					self.anim.buttons_up('liberty', self.liberty)
				

				# boutons
				elif bY <= ev.pos[1] <= bY + bH:
				
					# stop and play
					if bX <= ev.pos[0] <= bX + bW:
						if self.stop:	
							self.stop = 0
						else:
							self.stop = 1
						self.anim.buttons_up('play', self.stop)
						
					# H
					elif bX+bW+DECAL <= ev.pos[0] <=  bX+bW+DECAL+32:
						if self.H[2] == 0:
							self.anim.buttons_up('field', (cX, cY))
							self.H[2] = 0.0000001 
						else:
							self.anim.buttons_up('field', None)
							self.H = [0, 0, 0]
							# recalcule l'energie totale
							self.init_E_M()


					
				# initialisation de la grille		
				elif iY <= ev.pos[1] <= iY + iH:
					if iX <= ev.pos[0] <= iX + iW:
						self.init_grid('hot')
						self.anim.init_grid(self.spin)
					elif iX + 25 <= ev.pos[0] <= iX + iW + 10 * R:
						self.init_grid('alea')
						self.anim.init_grid(self.spin)
					elif iX + 50 <= ev.pos[0] <= iX + iW + 20 * R:
						self.init_grid('cold')
						self.anim.init_grid(self.spin)
					
						
	def monte_carlo (self, times = -STEP):
		" itere l'algo de monte carlo "

		while times:
			
			# on reactualise l'affichage tous les STEP iterations
			if (times % STEP == 0) and self.display :
			
				# on intercepte les evenements
			 	self.event_handler (event.get())
				 
				# on stoppe l'animation?
				if self.stop :
					time.delay(400)
					continue
				
				# raffraichi les graphes
				max_E = 2 + max(abs(self.H[0]), abs(self.H[1]))
				self.anim.graph_E.trace(self.anim.screen, self.E_total / float(NB_SPIN), max_E * 1.2)
	 			self.anim.graph_M.trace(self.anim.screen, self.M_total / float(NB_SPIN), 1.2)
				
				# affiche le temps de Monte-Carlo
				self.anim.buttons_up('tmc', -int(times / NB_SPIN))
				
				# raffraichi l'animation
				self.anim.refresh()
				 
				 
			# choix aleatoire d'un spin sur la grille
			i, j = randint(0, GRID[0] - 1), randint(0, GRID[1] - 1)
			
			# on choisit un etat aleatoire
			etat = self.get_new_spin(self.spin[i][j])
			
			# calcule la variation d'energie pour cet etat 
			#	en fixant les conditions aux limites periodiques
			dE = 0
			
			for n in xrange(3) :

				# on calcule la projection des spins voisins sur chaque axe (avec effets de bords)
				if self.bord :
					s = 0
					if i != 0 : s += self.spin[i-1][j][n]
					if i != GRID[0]-1 : s += self.spin[i+1][j][n]
					if j != 0 : s += self.spin[i][j-1][n]
					if j != GRID[1]-1 : s += self.spin[i][j+1][n]

				# sans effets de bords
				else :
					s = (self.spin[(i-1) % GRID[0]][j][n] + 
					self.spin[(i+1) % GRID[0]][j][n] + 
					self.spin[i][(j-1) % GRID[1]][n] + 
					self.spin[i][(j+1) % GRID[1]][n])
						
				# calcule la variation d'energie entre l'ancien etat et le nouvel etat
				dE += (etat[n] - self.spin[i][j][n]) * (self.J * s - self.H[n])
									 

			# si l'energie est abaissee on lui affecte etat (selon la methode de metropolis)
			if dE <= 0 or random() < exp(-dE / self.temp):
				 
				 # calcule l'energie totale
				 self.E_total += dE
				 
				 # calcule la magnetisation totale (somme des projections sur chaque axe)
				 self.M_total += etat[0] + etat[1] + etat[2] - self.spin[i][j][0] - self.spin[i][j][1] - self.spin[i][j][2]
				 
				 # affecte le changement d'etat au spin[i,j]
				 self.spin[i][j] = etat
				 
				 # on affiche le spin qu'on vient de modifier
				 if self.display : self.anim.change_spin (i, j, self.spin[i][j])
				 				 
			# on decremente
			times -= 1


if __name__=='__main__':
	

	#parametres de la classe IsingModel
	#IsingModel(style = {cold,hot,alea}, display = {0,1}, nb_etat = {0,2}, liberty = {1,2,3}, H = [x,y,0], J=-1, grille = [p,q])

	ising = IsingModel('hot', 1, 2, 1, [0,0,0], -1, [70,70])

	ising.monte_carlo()
